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ABSTRACT 

In this study an accurate finite-difference scheme has 
been utilized to investigate oscillatory, laminar and incompress- 
ible flow between two-parallel-plates and in circular tubes. The 
two-parallel-plates simulate the regenerator of a free-piston 
Stirling engine (foil type regenerator) and the channel wall was 
included in the analysis (conjugate heat transfer problem). 
The circular tubes simulates the cooler and heater of the engine 
with an isothermal wall. The study conducted covered a wide 
range for the maximum Reynolds number (from 75 to 60,000) , 
Valensi number (from 2.5 to 700), and relative amplitude of 
fluid displacement (0.714 and 1.34). 

The computational results indicate a complex nature of 
the heat flux distribution with time and axial location in the 
channel. At the channel mid-plane we observed two thermal 
cycles (out of phase with the flow) per each flow cycle. At this 
axial location the wall heat flux mean value, amplitude and 
phase shift with the flow are dependent upon the maximum 
Reynolds number, Valensi number and relative amplitude of 
fluid displacement. At other axial locations, the wall heat flux 
distribution is more complex. 

INTRODUCTION 

Recently, a great deal of work has been devoted to study 
unsteady fluid flow and heat transfer both computationally and 
experimentally. In this paper attention is given to one 
particular class of unsteady flows namely oscillatory flow with 
zero mean input velocity. This is in contrast to other unsteady 
flows that have non-zero mean flow and known as "pulsatile 
flow”. 

One of the applications of oscillatory flow analyses is in 
the free-piston Stirling engine heat exchangers. In these 
engines the working fluid oscillates from the compression space 
to the expansion space going through the heat exchangers 
during the process. In this cyclic process the working fluid 
absorbs and releases energy within the various heat exchangers 
resulting in a net work output to the power piston. The design 

( 


codes modeling these type of engines currently use steady state 
correlations for estimating the friction and heat transfer losses. 

In an oscillatory flow, the cyclic motion results in the 
velocity and temperature profiles being significantly different 
from those obtained from unidirectional or steady flows [1-3], 
which significantly effects the friction factor and heat transfer 
coefficient. Kurzweg [l] analytically solved the conjugate heat 
transfer problem for flow within two-parallel-plates channel, at 
the mid-plane of the channel (i.e. fully developed conditions), 
Kurzweg found that significant increase in the axial transport 
of heat can be achieved by flow oscillations without a signifi- 
cant mass transport between the hot and cold reservoirs. 
Similar analysis for circular tubes can be found in [4], 

In this paper attention has been given to the engine heat 
exchangers namely: the I) Regenerator, 2) Cooler. 3) Heater. 
The regenerator numerically modeled is the foil type which can 
be represented by a two-parallel-plates geometry (Ibrahim et 
al.[3]). Also, the cooler and heater are modeled as a single 
circular tube with constant wall temperature. 

Figure 1 shows the envelope in which different Stirling 
engines operate, plotted in terms of the Valensi number ( Va) 
and maximum Reynolds number ( and the cases 
investigated in this paper are also marked. The different 
criteria for transition from laminar to turbulent flow are shown 
by different straight lines [5,6]. For low Re^, (left side of the 
lines) the flow is laminar throughout the cycle; while for high Re m 
(right side of these lines) the flow is a combination of lami- 
nar/transitional/turbulent over the cycle. As it is evident from 
the plot most of the Stirling engines operate in the transition 
or "fully turbulent" zone. Ahn and Ibrahim [7] tried to map 
the zones wherein quasi-steady turbulence models could be 
used for oscillating flow conditions, they used a High Reynolds 
Number turbulence model and found the predictions to be poor 
in the transition region. Recently an empirical transition 
model has been proposed for oscillatory flows [8], this model 
hopes to improve the predictions and is currently undergoing 
some testing. 
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Since computational turbulent or transitional flow 
conditions require experimental data for comparison (these 
data are not available today), the present study is restricted to 
laminar flow conditions. 
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ANALYSIS 

Assumptions 


Figure 2 shows the two geometries, with the inlet and 
boundary conditions examined in this paper. Figure 2a shows 
the conjugate problem where the channel wall is included in 
the computation domain. In this configuration the two-parallel- 
plates channel is used for code validation purpose as well as 
modelling the regenerator (foil type), of the Stirling engine 
(cases Rj & R2) employing the Cartesian coordinate system. 
Figure 2d shows the circular tube problem where the computa- 
tion domain is limited to the fluid only, the tube wall being 
isothermal. This circular tube geometry is used to model the 
cooler (cases C^, C2 & C3) or the heater (cases H, , H2 & H3) 
of the Stirling engine employing the cylindrical coordinate 
system. The following assumptions are made: (1) the flow is 
laminar incompressible with constant thermophysical proper- 
ties; (2) The inlet velocity profile is uniform but oscillating 
with time; (3) As for the heat transfer inlet and boundary 
conditions two cases are distinguished: i- (Cases R^ & R2) the 
fluid enters the channel with uniform temperature hot from one 
end and cold from the other end, maintaining uniform axial 
temperature gradient. The fluid exchanges heat with the solid 
wall which is subjected to the same uniform axial temperature 
gradient, ii- (Cases Cj, C 2 , C 3 , H^, H2 & H,) the fluid enters 
the channel with uniform temperature and is heated (or cooled) 
by heat transfer from (or to) constant temperature walls; (4) 
The axial viscous diffusion as well as heat conduction are 
negligible. 

Mathematical Model 


The system of partial differential equations in Cartesian 
and axisymmetric coordinate systems used to model the 
unsteady flow are as follows: 


Continuity Equation: 


du 

dx 


♦ 


0 


( 1 ) 


12=1 for the axisymmetric coordinate system and, 
12=0 for the Cartesian coordinate system with r=y. 


Here S u and S v denote the source terms in the momentum 
equations and are given as follows: 
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and P is the hydrodynamic pressure. 


Energy Equation: 
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5 r is the source term in the energy equation and is assumed to 
be zero or negligible in the paper. 


Boundary Conditions: 

The boundary conditions used to solve the above equations 
are: 

a) At the walls: 

u=v*OanJ{T*T wall OR q" - q ) (*) 


b) At the axis of symmetry: 

— = v = 0 . q" - 0 (8) 

dr 

c) Inlet plane: 

. v fc - 0 , r* - ir^ (9) 

d) Outlet plane: Since the values of the dependent 
. variables are not known a priori at these planes the gradients 

of the dependent variables in a direction normal to the outlet 
plane are neglected. Such a situation is valid if the outlet plane 
is far from the entrance or any recirculating activities. 


Momentum Equations: 

The conservation of momentum yields two scalar equations in 
the x and r directions respectively. 


- E - E . o 

dx dx dx 
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x-Momentum : 



Numerical Method 


The governing equations are solved numerically by a 
conservative finite volume (FV) method utilizing a modified 
version of the computer code CAST, developed by Peric and 
Scheureur [9]. The code has been extended to include time 
dependent boundary conditions (Oscillatory flows) and solve 
for conjugate heat transfer type problems. The system of linear 
algebraic equations for each conservative equations are solved 
by the efficient strongly implicit procedure (SIP) (10J. The 
nonlinear iterations are done by using the SIMPLE algorithm 
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A 52x52 grid was used for all the cases investigated with 
grid density being high close to the walls and sparse away from 
it in the transverse direction, and a uniform grid was used in 
the longitudinal direction (axial). Convergence criteria was set 
at 0.1 % of the global residual norms for every dependent 
variable. For each case 120 time steps per cycle were used and 
at least 4 oscillation cycles were needed to reach quasi steady 
oscillating flow conditions. A typical run involved approxi- 
mately 1000 seconds of CPU time per cycle on a Cray YMP/24 
(sn 1040) supercomputer. 

RESULTS AND DISCUSSION 
CODE VALIDATION 

Kurzweg [lj analytically solved the N-S equations 
governing the oscillatory flow between two-parallel-plates 
channel. The geometry investigated by him is shown in 
Fig. 2a. The flow is set to motion inside the channels by a 
sinusoidally varying pressure gradient, and a constant linear 
temperature gradient is maintained along the channel through- 
out the cycle. Kurzweg assumed the channels to be long such 
that the fluid flow is fully developed. Under this assumption 
the momentum equations simplify considerably and are 
tractable to analytical techniques. 

In contrast to the analytical analysis, the numerical one 
has a two dimensional domain and finite; the flow is estab- 
lished by a sinusoidally varying velocity input 
[ u im 8« rt Sin(ut) ], as it is, relatively, more difficult to apply 
pressure boundary conditions numerically. Figure 3 shows 
comparison between analytical and numerical velocity profiles 
at different velocity phase angles for the conjugate problem, 
case R 2 : Rc laiu = 12000 and Va =400 (see Table 1). The 
symbols are used for the analytical solution [l], while the 
dotted lines are used for the present work; the agreement 
between the two is excellent. Here the velocity profile shows 
the presence of a small Stokes layer (high Valensi number case) 
with the flow field is almost uniform in the channel core. At 
some instants the velocity profile exhibits a flow reversal near 
the wall. 

Figure 4 shows comparison between analytical and 
numerical instantaneous temperature profiles ( T - T x ) at 
different velocity phase angles for the co 2 \jugate problem, case 
R 2 : Rc m . T = 12000 and Va =400. Here the temperature profile 


shows a steep gradient near the wall and almost flat in the 
channel core with a possible maximum near the wall at some 
velocity phase angles. It should be noted that at some velocity 
phase angles (e.g. 90) the temperature in the core is almost 
similar to that at wall. Using conventional heat transfer 
relationship this will mean that the heat flux at the wall should 
be zero. However, the results indicate the presence of a 
temperature gradient at the wall and accordingly a heat 
transfer at the wall. As can be seen from the plots the agree- 
ment between the analytical and numerical results is excellent. 

Cooler and Heater Cases 

Table 1 lists different cases investigated which cover a 
wide range of Va and A r . These cases as well as others, 

were discussed in details in Kannapareddy [13]. In this paper, 
for the space limitations, we will present case H 2 only . 

Figure 5 shows the temperature contours at different 
velocity phase angles for the circular pipe problem and case H 2 

Re JMX = 16500, Va- 88, L/ D b = 70 and A= 1.34. The 
contours are shown only for half of the cycle (0° to 180°) with 
cold fluid entering the tube (620 K) while keeping the wall at 
a hotter temperature (650 K). The entering fluid cold front 
advances into the channel during the acceleration portion of 
the cycle; this front continue advancing during the deceleration 
portion of the cycle but at a lower rate. Upon examining this 
case together with other cases we found that the velocity is out 
phase with the temperature; the degree of out phase is greater 
as the Re tMT (or Va) is higher. In contrast, as (or Va) is 
very low (e.g. case R*) the velocity and temperature are in 
phase (not shown in this paper). 

Figure 6 shows the normalized values for Nusselt 
number, wall heat flux and temperature difference at the 
channel axial mid-plane, for case H 2 : 16500 and Ka= 

88 . The temperature difference is obtained as the absolute 
difference between the wall and the section average fluid 
temperature on Fig. 6a and the bulk fluid temperature on 
Fig. 6b. On both figures the wall heat flux is the same, while 
the Nusselt number is different because of the temperature 
difference employed. 

Using the section average temperature, the thiee 
quantities described above (Nusselt number, waill heat flux and 
^temperature difference) are in phase with each other (Fig. 6&). 
On the other h.and, using the bulk temperature [12], the 
temperature difference (and accordingly the Nusselt number) 
is out of phase with the heat flux (Fig. 6b). Moreover, the 
temperature difference passes through zero and accordingly, the 
Nusselt number shoots to infinity. It should be noted that the 
two plots show symmetry in time i.e. the thermal cycle is 
repeated twice for one flow cycle, this is due to the inflow 
symmetry and the plots are made at the channel mid-plane. At 
other planes this symmetry with time is lost as expected since 
the flow enters the channel at a uniform temperature. 



Several other cases have been plotted (not shown in the 
paper) for other Va and A f . These cases were plotted at 

different values of x/D^. The plots reveal the complex nature 
of the heat flux distribution with time and axial location in the 
channel. Efforts are underway to Fourier analyze these results 
and try to obtain a heat transfer correlation of some practical 
use. 

CONCLUDING REMARKS 

In this paper, we computationally examined the fluid 
flow and heat transfer in three different components of the 
Stirling engine namely regenerator, cooler and heater. Some of 
the cases examined are summarized in Table 1 and they 
represent the operating conditions of the NASA Space Power 
Demonstrator Engine. 

Cases R^ and R£ resemble the regenerator and have 
been modeled using the conjugate heat transfer problem with 
a two-parallel-plate channel. Cases C^, C 2 & C 3 as well as Hp 
& Hj resemble the cooler and heater respectively and have 
been modeled using the circular pipe with an isothermal wall. 

The computational results revealed the following: 

1- For low Re ^ and (or) Va , the fluid flow and heat transfer 
for the regenerator (case Rj) are quasi-steady i.e. the velocity 
and temperature profiles are parabolic. For high Re ^ and (or) Va 
(case R 2 ) flow reversal takes place near the wall at some parts 
of the cycle; also almost flat velocity and temperature 
profiles are obtained in the core of the channel. These cases 
were not only used for studying the foil type regenerator, but 
also to validate the code by comparing the computational 
results with analytical solution at similar operating conditions 
(the comparison showed excellent agreement). 

2- The heat transfer results for the cooler and heater cases 
show that the heat transfer (in the channel mid-plane) goes 
through two cycles per each flow cycle. Also, the temperature 
profile is out of phase with the velocity profile particularly at 
the high Re m and (or) Va. 

3- The usual definition of the heat transfer coefficient in the 
case of oscillatory flows is ambiguous and limited due to the 
way the temperature difference is determined. The common 
practice has been to use the mixing cup or bulk temperature as 
the reference for the temperature difference as in the case of 
unidirectional or steady flows. But the bulk temperature 
definition (velocity weighted temperature across a cross 
section) breaks down in oscillatory flows due to flow reversal 
close to the wall for parts of the cycle at high Va and is 
especially true for laminar flows. A way of resolving this is to 
use the absolute value of the velocity in the definition of the 
bulk temperature and then evaluate the heat transfer coeffi- 
cient. In this study both the section average temperature and 
bulk temperature were explored. 


4- Using the section average temperature, the three quantities, 
namely, Nusselt number, wall heat flux and temperature 
difference are in phase with each other and out of phase with 
the flow. The wall heat flux is symmetric at 180 only in the 
mid-plane of the channel, other wise it has a different shape in 
the first half of the cycle compared with the second half. Also 
the wall heat maximizes at a lower velocity phase angle for 
larger A f values (maximum wall heat flux in the mid-plane of 
the channel takes place at 147° for case C~ and at 105° for case 

H 2 ). 

5- On the other hand, using the bulk temperature, the temper- 
ature difference (and accordingly the Nusselt number) is out of 
phase with the heat flux (Fig. 6b). Moreover, the temperature 
difference passes through zero and accordingly, the Nusselt 
number shoots to infinity. 

6- Several computational cases have been examined for other 
R e^ , Va and A r (not shown in the paper). The wall heat flux 
plots at different values of x/D^ reveal the complex nature of 
the heat flux distribution with time and axial location in the 
channel. Efforts are underway to Fourier analyze these results 
and try to obtain a heat transfer correlation of some practical 
use. 
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NOMENCLATURE 


a 

Half channel width in the conjugate 
problem. 

a. 

Relative amplitude of fluid motion: 
* (2X 1Blx ) / L 
1 1 te^ 

2 '(I/D*)' Va 

b 

Half solid thickness in the conjugate 
problem. 

c. 

Specific heat of the fluid or solid. 

D h 

Hydraulic diameter of the tube or channel. 

k 

Thermal conductivity. 

L 

Length of the channel or tube. 

Nu 

Nusselts number. 

P 

Hydrodynamic pressure in the momentum 
equations. 

Pr 

Prandtl number of the fluid. 

4" 

Heat flux. 

r 

Radial coordinate distance and direction 
for Axisymmetric cases or, 
the radius of the tube. 



K H 


t 

T 

T. 

T> 


Maximum Reynolds number: 

= <P- «„,•£*) / P 

Time [sj. 

Temperature in K, 

T(x,yj) or T{xf,t ) . 

Section averaged temperature of the fluid, 

W) ~\fTrdr]j\frdr] 

Bulk temperature of the fluid weighted by 
the absolute velocity, 

T b {xj) -[/ |«(r)| Trdr] / [f|«(r)| rdr] , 
Time averaged temperature, 

T x {xj) “[T^+yx] . 


u 



V 


Va 


y 

Greek : 

u 

Y 

v 

P 
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Axial velocity in the x coordinate 
direction , u (xyj) or u {x,rj) . 

Axial velocity at the inieL . 

Maximum input velocity in a cycle. 

Normal velocity of the fluid in the y 
coordinate direction. 

Valensi Number: 

Va = (p.u.D* 2 ) / (4.p). 

Axial coordinate distance and direction. 
Maximum amplitude of fluid displacement. 

Normal coordinate distance and direction 
for Cartesian coordinate system. 

Angular frequency in rad/sec. 

Dynamic viscosity of the fluid. 

* -(077 dx ) ■* / L ,time averaged 
constant axial temperature-gradient. 
Kinematic viscosity of the fluid. 

Density of the fluid. 
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Figure 1 Envelope in which different Stirling Engines operate 
together with; 

i) Criterion for transition from laminar to turbulent flow, 

ii) Different cases studied in the present work. 
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Figure 3 Comparison between analytical and numerical velocity 
profiles at different velocity phase angles for the conjugate problem 
for Case R 2 12000 and Var 400. 

(Symbol : Analytical- Kurzwcg [1] 

Dotted line : Numerical- Present work) 
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Figure 2 Different geometries, inlet and boundary conditions 
examined; 

a) Conjugate problem, flow between two parallel plates, 

b) Flow inside circular tube. 
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Figure 4 Comparison between analytical and numerical temperature 
profiles at different velocity phase angles for the conjugate problem 
for Case R 2 R e^ - 12000 and Var 400. 

(Symbol : Analytical- Kurzweg [1] 

Dotted line : Numerical- Present work) 
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Figure 5 Temperature contours at different velocity phase angles 
for the circular pipe problem, Case H 2 16500 and Va m 88. 

[Temperature in K : 

A:649.50, B:640.00, C:630.00, D:625.00, E62250 

F:620.90, G:620.30, H:620.08, 1:620.01, J^2G£3qi 


Vdoaty Phase Angle 

■ = Nu/Nu,max 
d= Q/Q,max 
o = (Tb— Tw)/(Tb-Tw)max 

Figure 6 Normalized values for Nusselt number, wall heat flux and 
temperature difference at the channel axial mid-plane, for 
Case H 2 - 16500 and Va - 88, using: 

a) Section average fluid temperature (7 fl )* 

b) Bulk fluid temperature (7*). 









